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The lattice formulation provides a way to regularize, define and compute the Path In- 
tegral in a Quantum Field Theory. In this paper we review the theoretical foundations 
and the most basic algorithms required to implement a typical lattice computation, in- 
cluding the Metropolis, the Gibbs sampling, the Minimal Residual, and the Stabilized 
Biconjugate inverters. The main emphasis is on gauge theories with fermions such as 
QCD. We also provide examples of typical results from lattice QCD computations for 
quantities of phenomenological interest. 
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1. Introduction 

In this paper we present a brief review of Lattice Quantum Field Theory (LQFT) , 
a way to formulate a Quantum Field Theory (QFT) in algorithmic terms a . 

Most of this work is based on lectures given at Fermilab by the author in 2001. 

QFTs are the application of quantum mechanics to fields. They form a very 
general class of mathematical models that reduces to quantum mechanics in the non- 
relativistic limit (speed of light — > oo), to relativistic mechanics in the decoherence 
limit (Plank constant — > 0) and to classical physics when both limits are taken. 

Any QFT states that b : 

• Each type of elementary particle "A" is associated with a field <fi(x) so 
that |0(a;)| 2 represents the probability of the event "particle A is at the 
space-time point x — (a;n,x)"; 

• There exists a functional of the field S[<j), ...], called action, which describes 
the dynamics of the system and any correlation amplitude between multiple 
time ordered events A at x, A at x' , A at x" , etc. and can be computed as 



The square of a correlation amplitude is a regular real correlation function. 



a For an introduction to QFT see l > 2 and for LQFT see 3> 4 .5>6. 

b From now on we assume units in which both the Plank constant and the speed of light are 1. 
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The field components <p(x) can be scalars, complex, spinors, vectors, tensors, etc., 
depending on the properties of the particle A. x', x', x" are points in the space-time. 
The symbol J[d4>] indicates an integral over a Hilbert space and is known as Path 
Integral (PI) . Each (f> in the Hilbert space is referred to as a history or path or field 
configuration. 

The symbol <f> is used here as a template for a generic field and, for now, one 
may think of it as a real scalar field. Later <f) will be replaced by the symbol U when 
it represents a gauge field and by the symbol tjj when it represents a spinor, for 
example, a quark. 

Translational invariance combined with the finite speed of light implies that S 
can be written as the integral over the s-dimensional space-time of a local function 
of <fi, the Lagrangian density C. 

S[4>,...] = J £(Hx),d^(x)...)d s x (2) 

The dots indicate additional fields and/or parameters of the theory such as 
particle masses and coupling constants. The Lagrangian is usually a function of 
4>{x), it's gradient and higher derivatives. 

The naive assumption that the Hilbert space in the PI is the space of all possible 
distributions leads to the problem of divergences of QFTs. In order to understand 
and cure this problem one has to properly define this space and provide an algo- 
rithmic way to evaluate the integral. That is the main purpose of these notes. 

For some theories, such as quantum electrodynamics (QED), it is possible to 
perform a functional expansion of the integrand of eq. (1) around a minimum and 
integrate exactly the individual terms. These terms are the Feynman diagrams. The 
result is an asymptotic series, the Dyson series, that can be used to approximate 
the PI. While this is a well defined algorithm it is very difficult to implement, it 
only works in some cases, and it cannot be improved to arbitrary accuracy because 
of the asymptotic nature of the Dyson series. The divergences that appear in this 
case can be subtracted by regulating the theory. For example, one can dump the 
Fourier modes of the fields with frequency above some cut-off p. 

The perturbative expansion plays a role analogous to the Taylor expansion and 
it does not provide a general definition of the PI. Moreover, for some theories, the 
Dyson series may diverge and a different approach is required in order to define and 
compute eq. (1). 

For quantmm chromodynamics 7,8 (QCD) and other strongly coupled theories, 
the lattice formulation has been one of the most successful. We will refer to the 
latter as lattice quantum chromodynamics (LQCD). 

2. Overview 
2.1. Formulation 

The formulation of LQFT consists of three stages 3 ' 4,5 : 
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Step 1: Discretization. The space-time is approximated with a finite hyper- 
cubic mesh such that <p(x) is only defined on the lattice sites x. Now the Hilbert 
space that constitutes the integration domain is well defined. The hypercubic mesh 
is characterized by the number of dimensions, s, the lattice spacing, a, the overall 
size, L = Ka, and the boundary conditions. Here is an example of a 3D lattice. 



<j>(x) at site x ■ 



• • • • 



L 



a 



After discretization, for every finite a the symbol f[d(j>] becomes a well-defined 
multidimensional integral 

[d<f>] - / #o / d<h Jdfa... (3) 

where the i-th degree of freedom, fa is localized at some lattice site. 

Step 2: Wick rotation. The time coordinate xo is Wick rotated x — ► ix thus 
turning eq. (1) into 

(0| \0) E = / [d<f>]<f>(x)<f>(x')cj>(x")...e- s °W (4) 



where Se is now the Euclidean action. If the exponent term is real, and this is 
true for most physical systems of practical interest, then eq. (4) reads as a weighted 
average of <j>{x)4>(x' )4>{x" ') . . . with an exponential weight factor equals to exp(— Se)- 

Eq. (4) looks like a correlation in statistical mechanics and the integral can be 
computed using standard numerical technques. 

Step 3: Monte Carlo Computation Eq. (4) is approximated by a finite sum 

, fe<JV 

(0| <p(x)<p(x')<p(x")... |0> B <P [k] (x)<t> [k] (x')<j> [k] (x")... (5) 

fc=0 

Each field configuration <p^ is generated by a random sampling from a prob- 
ability distribution proportional to exp(— >Se[<^]). As more terms in the sum are 
considered, the right hand size of eq. (5) approaches the left hand side. The nu- 
merical error in this numerical approximation can be controlled and it is usually 
proportional to N~? where TV is the number of generated configurations. 

This lattice approach to the PI can be improved to any arbtrary precision by 
reducing the discretization step (a — > 0), by considering a larger portion of space- 
time (L — ► oo), and by generating more configurations (N — ► oo). The extrapolation 
of a to zero is referred to as the continuum limit. 
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1 Algorithm: Metropolis for MCMC 

2 Input: 4>, Euclidean action S 

3 Output: </>' 

4 

5 generate <f>' at random uniformly in the Hilbert space 

e generate a random number z 

7 if z > exp(S(<f>) -5(0',...)) then 
s replace <f>' with 4> 

9 return <t>' 



Fig. 1. Metropolis is the simplest algorithm to generate a Markov Chain Monte Carlo (MCMC). 
It performs a global change of the input configuration and then an accept-reject step of that 
change. The latter step ensures the reversibility condition and the crgodicity of the chain. Notice 
that every time an update is false the output configuration is the same as the input. 



Notice that the set of Monte Carlo configurations which appear in the sum 
eq. (5) does not depend on the particular correlation amplitude one is computing 
and therefore it can be reused for different computations as long as the system is 
the same (i.e. it is described by the same action). For example, in typical LQCD 
computations, a large numbers of configurations are generated, stored, and reused 
for different computations. 

Because of the Wick rotation, any correlation computed using the above tech- 
nique is defined in the Euclidean time, not in the Minkowsky time. Some observables 
are unaffected by the Wick rotation and can be reliably computed, while others re- 
quire an analytical continuation back to Minkowsky time. Some phase information 
may be lost when combining this analytical continuation with the finite precision of 
the Monte Carlo method and only those observables that do not require analytical 
contination back to the Minkowsky time can be reliably extracted from the lattice 9 . 
Luckily these include the low energy spectrum and the absolute value of matrix ele- 
ments of operators between on-shell states. In LQCD typical computations include 
the masses of bound states such as glueballs, mesons, baryons, and matrix elements 
between these states. 

2.2. Algorithms 

Any Monte Carlo computation can be decomposed into three elementary steps: 
Step 3.1: Markov Chain Monte Carlo (MCMC). 

The MCMC is a method to generate the field configurations </> by random sam- 
pling from a known probability distribution which, in this case, is proportional to 
exp(-S E [4>})- 



5 



3 



2 



Algorithm: Gibbs sampling for MCMC 
Input: 4>, euclidean gauge action S 
Output: 4>' 
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5 for all lattice sites x 

e store the field variable </>(x) in 4> old 

7 replace the value of 4>(x) with a random one (uniform in domain) 

8 generate a random number z 

9 if z > exp(— [change in action]) then 

10 replace <j>{x) with 4> old 

11 = <j> 

12 return 6' 



Fig. 2. Gibbs sampling is another algorithm to generate the MCMC. It loops over all degrees of 
freedom of the system and, for each of them, it performs a local change and an accept-reject of 
the change. If the action is local, Gibbs sampling and algorithms derived from it are more efficient 
than the Metropolis. 



The idea behind the Markov Chain is that of bulding a randomized iterative algo- 
rithm M. so that the transition probability Pm{$\<$) of going from a configuration 
4> to a configuration <fi' satisfies the following reversibility condition 



Regardless of the starting configuration </>[ ] , the succession in k is ergodic with 
the desired stationary distribution 10 exp(— 5^(0, ...)). The simplest MCMC algo- 
rithm that satisfies the reversibity condition, eq. (6), is the Metropolis 11 algorithm 
shown in fig. 1. This algorithm makes the next configuration in the chain by either 
picking a copy of the preceding one, or a totally new configuration generated with 
uniform probability in the configurations' domain. This choice is implemented as 
an accept-reject step that depends on the action and guarantees that the transition 
probability satisfies the reversibility condition. 

One algorithm derived from the Metropolis is the Gibbs sampling, fig. 2. It uses 
an accept-reject similar to the Metropolis algorithms but differs because only one 
field variable is updated at one time. If the action is local, and this is usually true 
for most physical systems, the accept-reject condition is also local, and the overall 
algorithm is more efficient than the Metropolis in sampling the space. 

There are many other algorithms that can be used for generating the MCMC 
configurations and most of them are derived from either the Metropolis or the Gibbs 



<t>[k\ — M <t>[k+i\ 



(6) 
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sampling. 

Step 3.2: Measure the operator. This algorithm measures the desired cor- 
relation, the operator O = <f)(x)<f)(x')(fi(x")..., on each MCMC configuration <f>^, 



4>[k\ ■* o 



<t>[k](x)<t>[k](x')<t>[k](x")~- 



This step is non-trivial because the operator O may depend on fields that do not 
appear in the PI, for example fermions propagating in a background gauge field. 
If this is true, the algorithm O requires the computation of fermion propagators. 
We'll provide more details in a later section. 

Step 3.3: Average and Analysis. 

The final result is computed by averaging the measurements of each configura- 
tion. This average is accompanied by a statistical analysis to determine the error in 
the result. 



<A[i] (^^(y^ii]^")--- 

4>[2] (x)4>[2] (x')<l>[2](x")- 
<t>[3] O)0[3] (x')<t>[3](x")--- 




(O|0(x)0(xW)...|O> 
± statistical error 



Estimating the stastical error is of crucial importance since it must used as a stop- 
ping condition for the MCMC. A naive estimate of the error is given by o j\f~N where 
a is the standard deviation of the individual measurements Oua. However, this es- 
timate fails when the individual measurements are not Gaussian distributed. The 
standard algorithm to estimate the statistical error in an average without making 
any assumption about the underlying distribution is the Bootstrap algorithm 12 . 

Every lattice computation is a combination of the above elementary steps as 
shown below. 



Se, 4>[o] ■* 

physical ■* 
irrelevant -* 
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■ 4>\2] — <\ O 

■ 0[3] ■ 



o 



continue MCMC .. 



A 



(Q\4>{x)4>{x')...\0) 
± statistical error 



Because of the way the configurations (p^ are generated, each of them retains 
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some memory of the preceding configurations in the chain. This correlation dies out 
with the distance in MCMC steps between configurations. If the evaluation of the 
operator O on each individual configuration is more expensive than the elementary 
MCMC step M. , it may be wise to skip a number of configurations and evaluate the 
operator only on a subset of the total number of MCMC configurations. By skipping 
c configurations at each step, for c large enough, <j>^ and 4>[k+ c ] win be sufficiently 
decorrelated and the procedure provides a better sampling of the integration space 
while saving CPU time. An empirical choice is making c larger than the maximum 
distance between lattice sites in lattice units, c > sL/a. In fact, if one assumes that 
for a local action, at each elementary MCMC step information propagates only from 
one lattice site to the next, after c steps, information should propagate all around 
the lattice, thus removing the memory of the preceding configuration. 

From now on, when referring to a MCMC step, we will indicate the repetition 
of multiple elementary steps in order to achieve a sufficient decorrelation between 
effective consecutive configurations. 

2.3. Input of The Algorithms 

Any lattice computation takes the following input: 

• The action Se- The action determines the dynamics of the system. The 
same system may be represented by multiple actions that differ from each 
other because of irrelevant operators. These are high dimensional operators 
that can be added to the the Lagrangian and whose contribution to the 
action and to the correlation amplitudes vanishes in the continuum limit. 
Nevertheless, their contribution affects the rate of convergence of correla- 
tions to the continuum limit and can be relevant from a numerical point of 
view. The art of finding these operators and determining their coefficients 
is called improvement and is based on the work of Symanzik 13 . 

• The initial field configuration, <f>^ , from which the Markov Chain is started. 
The result of the computation should be independent from this choice be- 
cause the MCMC loses memory of this initial configuration. One way to 
verify that this is true is by repeating the computation using different </>[o] . 

• The lattice volume, represented by the parameter L. Ideally, one would like 
to perform computations close to the limit L — > oo. In practice, this is 
not possible. A finite L acts as an infrared cut-off which results in system- 
atic errors known as finite size effects. These effects can be quantified by 
repeating the computation at different values of L. For large L boundary 
conditions become irrelevant but, for finite L, they do affect the compu- 
tation. The usual choice consists of periodic boundary conditions in the 
hypercubic lattice topology: 



Vx, 4>{x + L) = J <f>(x) 



(7) 
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For fermions often antipcriodic boundary conditions are adopted 

Vx, tp{x + L) d = f -ij>(x) (8) 

For typical LQCD computations L ~ 32a ~ 2/m and often L is chosen 
larger in the time direction than in the space direction. Empirical results 
indicate that for L > 5/m w the corresponding systematic error is less than 
1%. 

• The number of MCMC steps N. This number is limited by the total com- 
puting time available. For different operators, O, the Monte Carlo integra- 
tion converges at different rates with N. For typical LQCD computations 
N ~ 1000. 

• Physical parameters of the Lagrangian, such as masses of elementary par- 
ticles, to, and coupling constants, g. Particles with masses m < 1/L and 
m > 1/a have propagators that, in Euclidean space, exhibit a correlation 
length larger then the lattice size L and smaller than the lattice spacing a, 
respectively. Elementary particles with these masses cannot be put on the 
lattice in a naive way. The standard approach for light fermions consists 
of performing the computation with non physical values for the masses (in 
the allowed region 1/L < m < 1/a) and then extrapolating the results to 
the physical values to — > to < 1/L. This extrapolation is called chiral ex- 
trapolation. For heavy fermions the extrapolation is just one of the possible 
approaches and different solutions are discussed in a later section. 

• Irrelevant parameters. One has the freedom to add irrelevent operators to 
the Lagrangian, i.e. operators whose contribution to the action vanishes in 
the continuum limit. Some choices for the coefficients of these operators are 
better than others because they improve the rate of convergence of the cor- 
relation amplitudes to the continuum limit. The values of these parameters 
do not affect the result of the computation but they do affect how fast one 
gets to the result within the required precision. 

Most notably, the value of a is not a direct input of any lattice computation. All 
the physical input parameters (masses and coupling constants) are bare parameters 
that, at constant physics, depend on the lattice spacing. Because of this implicit 
dependence, the choice of the coupling constant (which we'll refer to as g) is equiv- 
alent to setting the value of a. In general the relation betwen g and a is described 
by the Renormalizion Group Equation (RGE). 

Because the RGE is often only known perturbatively, one does not exactly know 
a priori the value of a in a lattice computation. Since every lattice quantity is 
computed in units of a, any error on a introduces an uncertainty in those quantities 
with dimension different from zero. For example the energy spectrum is in units of 
1/a. This problem is solved by computing ratios of masses and/or matrix elements 
that cancel any explicit dependence on a. An equivalent approach used in LQCD 
consists of measuring a by comparing the pion mass, m v , computed from the lattice 
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in lattice units with the physical pion mass and using this value to convert other 
quantities in physical units. We will see in the next section how this procedure is 
equivalent to choosing a renormalization condition. 

The RGE indicates that the contiuum limit a — > corresponds to a fixed point of 
the theory g — > g* (for example for asymptotically free 14 theories g* = 0) therefore 
the continuum limit is realized numerically by tuning g closer to the fixed point g* . 
Ideally the output of the computation should plateau as this limit is approached. 
In LQCD a typical example is the computation of the mass of the scalar over the 
vector meson for fixed bare quark masses. 

2.4. Regularization and The Continuum Limit 

On one hand the lattice provides a regularization scheme for the continuum PI. On 
the other hand the lattice formulation in the continuum limit provides a definition 
of the continuum PI. We wish to clarify the meaning of the concept of regularization 
and renormalization in the lattice paradigm. 
We start by making two observations 15 ' 16 : 

• One can never measure the value of a continuum field (p(x) at every point 
in space-time. One can only measure its integral over the test function that 
corresponds to a physical detector, which has a finite extension. Therefore, 
there are mathematical reasons to require that a continuum field be defined 
in the space of distributions. In particular, one assumes that a particle can 
be localized to any arbitrary precision and that the corresponding field 
configuration can be a delta function S(x). 

• One wants to model the short distance physics by introducing a Lagrangian 
density that contains only local (contact) interactions. Any non-trivial La- 
grangian contains terms that are products or powers of fields. 

These two observations are incompatible because the product or power of fields is 
not a well-defined quantity when a field is a delta function. The role of regularization 
is that of defining this product. 

There is a physical reason behind this problem: if one only knows the field via a 
finite size detector, with a spatial resolution limited to a, then any field fluctuation 
on a scale smaller than a should not be part of the model. This is why the model 
itself forces one to introduce some kind of cut-off a < a. The effect of modes with 
length scale smaller then a is encoded in the value of the bare coupling constants 
that one puts in the model. Any change in a implies a change in the modes that 
contribute to the bare coupling constants and, therefore, their values have to be 
scaled accordingly in order to describe the same physical system. 

The easiest way to regularize a distribution 6(x) is by replacing it with some 
localized function with finite support, for example 

S(x) -> 6 a (x) = i [6(x + a/2) - 9(x - a/2)} (9) 
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Fig. 3. [left]: Example of possible rcgularization for the delta function. The x axis is in units of 
a, the y axis is in units of 1/a. [center-right] : Example of a continuum field configuration <j> and its 
approximation with a linear combinations of regularized delta functions. 



This makes the product of distributions 5™(x) well defined. 

Let's consider now a ID scalar theory defined in the interval [0,L]. Its most 
general correlation amplitude, defined in terms of the PI, is 

(0\4>{xMx'Mx")...\0) = f J [d<t>]F[4>(x),g] (10) 

where F[<f>(x),g] is the integrand 

F[<j>(x),g] d ^ 4>{x)4>{x')(j>{x")...e- s ^^ (11) 

and g is a coupling constant that appears in the action. The specific form of the ac- 
tion, Se[4>i9] is unimportant but we assume that the action contains an interaction 
term of the form 

g J r(x)dx (12) 

with n > 2. This makes F[<f>,g] a non-trivial function of <fi. 

Lattice regularization is the way to regularize the integral (10) by approximating 
the fields with sums of regularized delta functions. 

K-l 

<f>{x) ~ ^ u (a, x) = f Y, <t>kta(x - ka) (17) 

k=0 

where 

-i r-ka+a 

<t>k d = - / Hx)dx (18) 

a Jka 

This is equivalent to discretizing the space-time on which the fields are defined (as 
shown in fig. 3). 

After discretization, for every finite K = L/a, 

J ( \d<j>]F[<j>(x),g] = f J d^ J d^...J d^_! F[0 latt (a, x), g] + 0(a) (19) 

S v ' 

K~L/a 

Here we used the upper index (a) , which identifies the regularized PI with lattice 
spacing set to a. 
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For a = 0.1 (K = 10): 

r(o) 



| [d<f>]F[<j>, g] = f J d<j> ... J d0 9 F[0, g] = 



For a = 0.05 (K = 20): 




+ ... + F[ 



J ° [&4>]F[^g] = f J d^,...y d0 19 i^,.g] = 





+ + 



For a = 0.025 (K = 40): 

r(o) 



y [dflFfag] = f J dfo-J d<j> 39 F[<j>,g] 





And at the continuum limit, a — > (if — > oo) 

F[0,5f] 



[cty]F[0 )fl ] =' lim 



K~L/a 




+ ... + F[ 



(13) 



+ ... 



(14) 



+ 



(15) 



+ ... 



(16) 



] + 



J <fi n (x)dx is finite f <p n (x)dx is finite f 4> n (x)dx~^oc 

Fig. 4. Example of lattice regularization of the PI 

Divergences associated with the limit a — > 0, K — > oo (at L = if a ^constant) 
are called ultraviolet, while those associated with the limit K, L — > oo (at a = 
L/if =constant) are called infrared. 

Fig. 4 shows, in a schematic way, how the path integral, eq. (16), can be approx- 
imated by finite multidimensional integrals, cqs. (13-15), and the fields are defined 
on a lattice. For each finite a, the "sum over the paths" (eqs. (13-15)) is well defined 
since all divergences that may appear can be absorbed in the normalization of the 
integration measure and, for each path <j>, the functional F[<j>, g] is finite. In the limit 
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a — > those configurations that correspond to a localized particle become more and 
more peaked and approach a Dirac S(x) function (eq. (16)). Since the integrand 
F[4>{x),g] is non-linear, the integrand diverges on those configurations cp(x) ~ S(x). 

In order to have a well defined limit a — > 0, if — > oo (at L = Ka =constant) one 
must require that the result of the regularized path integral is independent from a. 
The only way to do it is to making the field normalization and the coupling constant 
(the g of eq. (12)) dependent on a 

g^g R (a,A) (20) 

(the constant A must be introduced because, in general, a and g do not have the 
same dimensions). This makes the physics independent by the lattice scale a. 

One does this by choosing a particular correlation amplitude (identified by the 
functional integrand F[cj),g]) and imposing the contraint 



^0 (21) 



Eq. (21) is the RGE for a lattice regularized theory and it determines the behavior 
(the running) of gii{a). The appearance of A is called dimensional transmutation. 

A is the typical length scale of the physics being studied. This scale is in nature 
and there is no freedom to change it. For QCD, for example, it is best determination 
is from the LQCD scaling of the coupling constant 17 , A£f<f D ~ 259(l)(20)MeV 
(1/Aqcd * Ifm). 

Notice that this procedure of defining the limit a — > cannot be carried out for 
an arbitrary theory since there may be more sources of divergence than coupling 
constants. If this limit can be defined the theory is said to be renormalizable, if not, 
a must be kept finite and the theory should be considered as an effective theory. 

We distinguish between the bare parameters that appear in the regularized La- 
grangian (for a finite value of the cut-off, a) and the dressed or renormalized pa- 
rameters that are measured by actual experiments. If one takes the limit a — > 0, 
the bare parameters lose any physical meaning and one must carefully define the 
renormalized ones (one is said to choose a prescription). If one is happy with keep- 
ing the cut-off small but finite one is allowed to identify the renormalized and the 
bare parameters, because these can now be measured. This approach is known as 
the Kadanoff- Wilson approach to renormalization. 

In LQCD eq. (21) is realized numerically. One repeats the computation of the 
same quantity, for example the pion mass, m„, for different values of the coupling 
constant g 1 , g" , g'" , etc. thus obtaining a'fh w , a"fh w , a'"fh w , etc. where rh^ is the 
physical pion mass. By comparing the lattice results in lattice units with the physical 
pion mass one can obtain the value of a that corresponds to the input values of 
g. From the plot one determines g(a) and identifies the fix point g*. The same 



d 
da 



d^ / d<P 



d0 K _ 1 F[^(a;), 5fl (a,A)] 



K~L/a 
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Fig. 5. Different behavior of the integrand of b' n for high frequency modes (left) and low frequency 
modes (right) respectively. The x axis is in units of a. 

procedure can be carried on any physical quantity and it should be lead to the 
same running of g, up to corrections in a. 

Since we will never be able to probe the physical world at every length scale, 
every quantum field theory should be considered an effective theory. 



2.5. Lattice Regularization vs Momentum Cut-off 

A continuum field (f> can be expanded into Laplace components as 



<p{x) = Ke v - 



where 



n=0 



del 27m de} 1 

Pn = — ; K = 



4>{x)e 



c dx 



(22) 



(23) 



Similarly a lattice field, eq.(17) can also be expandend in Laplace components 
and we obtain 



where 



ilatt 



K-l 



(a, x) = J2 b'y 



71 = 



(ka - x)e- ipnX dx 



(24) 



(25) 



It becomes evident that for p n > 1/a the integrand oscillates quickly and the cor- 
responding integral, b' n , is small; while for p n < 1/a the integral is almost constant 
and approximately equal to e ~ lp " fca , therefore b' n ~ b n . The different behaviors of 
the integrand are shown in figure 5. This proves that eq.(22) can be written as 



4>(x) ~ latt (a, x) ~ <f> co (a, x)= f Y,8[-- Pn K 

n=0 \ a ' 



(26) 



The superscripts "latt" and "co" are used to identify the lattice and and the 
cut-off regularization schemes, respectively. 
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In other words the lattice cut-off a is equivalent to a momentum cut-off p max < 
1/a, an ultraviolet cut-off. Therefore the lattice and the momentum cut-off arc 
alternative but equivalent ways to regularize the PI. 

Note that 60 is the mean value of <j>(x) and p\ — 2ir/L represents the minimum 
energy /momentum mode that can propagate on a finite unidimensional volume of 
length L. This is a lower limit, an infrared cut-off. 



3. Pure Gauge Theories 
3.1. Action 

We consider a pure gauge theory and we restrict the gauge group Q to U(l) and/or 
SU{n). In order to be able to probe the gauge field we assume to have a single 
particle ip of infinite mass and unit charge that we can move on the lattice as we 
please. 

The Aharonov-Bohm experiment 18 teaches us that the gauge field is not directly 
measurable but if we move the test particle tp from one point 1 to a point x' along 
a path C, the test charge acquires a phase that can be measured via interference 
experiments 

il>(x r ) = exp (ig J A^x)dx^ $(x) (27) 

where A^(x) is the gauge field in the continuum space. On a lattice the shortest 
possible path is the link connecting two neighbor sites x and x+p, (+/t here indicates 
a positive vector in direction /i and length equal to the lattice spacing a) therefore 
the elementary lattice gauge degrees of freedom are not A^ but 



U(x,n) d = cxp (ig J A^(x)dxA ~ 1 + iagA^x + ^fi) + 0{a 2 ) 



(28) 



For later convenience we also define 

U(x, -n) d = f exp (igj^ A^dafj ~ 1 - iagA^x - + 0(a 2 ) (29) 

From the definition it follows that U(x, —fi) = W(x — fi, fj,). 

For the rest of this section the gauge links U(x, fi) will replace our generic tem- 
plate field <fi(x). The index \i is an integer that labels the vector component of the 
gauge field. 

Under a local gauge transformation i[)(x) — > V{x)ip{x) the field U transforms as 
U(x,n) -» V(x)U(x,v)V- 1 (x + fi) (30) 
Any path C on the lattice can be identified with a set of links U (x^ , ^j) such that 

exp J A' i (x)dx li ~ JJ?7(a;[ i ],/i[i]) (31) 
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The most local gauge invariant one can build is called the plaquette and it is a 
product of 4 links in a loop 

P(x, (j,, v) d = Tr [U(x, (j,)U(x + p,, v)U(x + /} + u, -fj,)U(x + v, -v)] (32) 
Here are some examples of plaqucttes: 



1* p 

P^l, 2) — £J.. .... 



The simplest lattice action Se one can engineer using the gauge field must 
therefore be linear in P(x,^,f), real, and invariant under the lattice rotational 
symmetry. Following the common notation this simplest action can be written as 

yauge dg zl Re P{X,H,U) (33) 

where n is the size of the gauge group (n = 1 for [7(1)). 
Substituting eqs. (28-32) in eq. (33) one obtains 

s9 au ge „ ^ ReY/x ^ [1 + lag A,(x + |A)][1 + iagA^x + \i>)] x 

[1 - iagA^(x - ^£)][1 - iagA u {x - ^z>)] 

■ 2 ■ 2 

■ 9 ■ 9 

[1 - ia^A,, - -^^[l - tag A, + —^-d^A v ] 



^£ x ^ 1 - ^(d„A v d^+g^A^f 



-c+^E^^f V" (34) 

where c is an overall irrelevant constant and eq. (34) was derived from eq. (33) via 
a Taylor expansion around x + \\i + jv. 

Notice that \F lll/ (x)F tlv (x) is the ordinary kinetic term for a continuum gauge 
field and a 4 ^ is / d 4 x in the contiuum limit. Following this analogy 

f> = ^ 05) 
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hence g is interpreted as the bare gauge coupling constant at the lattice scale a. For 
non-abelian gauge theories asymptotic freedom dictates that g goes to zero when a 
goes to zero. 

Eq. (33) is called Wilson gauge action 19 The Wilson gauge action can also be de- 
rived by direct discretization of the continuum gauge action by ignoring everything 
but the lowest order terms in a. 

For large (3 the trace of the average plaquette is close to n, perturbative effects 
dominate and this results in small quantum fluctuations, and small/slow changes in 
the configurations generated by the MCMC algorithm. For small fi the average pla- 
quette is close to 0, non-pertubative effects dominate which results in large/non-local 
quantum fluctuations, and relatively big changes in the configurations generated by 
the MCMC. 

3.2. Algorithms 

Because of the locality of the Wilson gauge action, for every link U(x,/j,) one 
can rewrite the action as 

ggauge = z£ ReTrU ( x ^ ^Rfc + ... ( 36 ) 

where the dots represent the sum over plaquettes that do not include U(x, /x) and 
R{x, /j.) = U(x + fi, v)U{x + fi + u, —fi)U(x + v, —v) + 

U(x + fi, — v)U{x + fi — z>, —n)U(x — i>, v) (37) 
are referred to as staples. A 3D projection is represented in the image below 



■ 

•• • 



This makes the Gibbs sampling algorithm more efficient than the Metropolis 
because it is possible to change a single link U(x,fi) — U old — > U(x,[i) = JJ new at 
one time without having to recompute the entire action. In fact, the accept-reject 
step only depends on the variation of the action given by the first term in eq. (36) 

6S gauge = zl ReTr yjjnew _ U<*d)) R fa ^)] ( 38 ) 

There are many algorithms that are similar to the Gibbs sampling but more 
efficient. One of the most common is the heatbath 20 . 
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Algorithm: Generate a random SU(n) matrix 
Input: n 
Output: A 

if n == 1 return exp(2ir i uniform^)) 
A = nx n identity matrix 
for i = to n - 2 
for j = i + 1 to n - 1 

a = 7r uniform () 

cp = 2 uniformQ 

cos(0) = 2 uniform() — 1 

sin(0) = a/1 -cos 2 (6>) 

u° = cos(a) 

u 1 = sin(a) sin(0) cos(</>) 
u 2 = sin(a) sin(0) sin(</)) 
u 3 = sin(a) cos(0) 
G = u° + ?i 1 cr 1 + u 2 cr 2 + u 3 er 3 
A' = i4 

for fc = to n - 1 

j^l ik A^ -\- G®^ A^ 

A' i k = G 10 A lk + G xl A^ k 
A = A' 
return A. 



Fig. 6. General algorithm for generating a random element in a U(l) and/or SU(n) gauge group 
with uniform distribution within the group. For SU(2) it uses the invertible map in 50(3) and 
for SU (n > 2) it generates the matrix as product of SU (2) subgroups. The function uniformQ is 
assumed to return a uniform random number in (0, 1). 



In order to make a MCMC step, whether global in the Metropolis or local in the 
Gibbs sampling and derived algorithms, one must be able to generate a new link at 
random with uniform probability in the gauge group Q. 

For Q = U(l) it is sufficient to generate a uniform random number u e [0, 1] 
and update 

U(x 7 ji) — > U'(x,[a) = e 27 ™ u (39) 

For Q = 5(7(2) one can use the map between SU(2) and SO(3) (the symmetry 
group of a 3-sphere), generate a uniform point on a 3-sphere (u°, u 1 , u 2 , u 3 ) and 
map it back into SU(2) via 

U(X, fJ,) -» U'{X, H) = U° + ?i 1 cr 1 + U 2 cr 2 + ?i 3 cr 3 (40) 

where Ui are the Pauli matrices. 
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Algorithm: Compute the static quark-antiquark potential 
Input: /3, size of gauge group n, number of MCMC steps N 
Output: V(r) 

Create local array V{r) and initialize it to zero 

for each lattice site x 
for each direction ^ 
set U(x, a) to a random element of the gauge group SU(n) 

for c = 1 to N 

replace U with the next configuration in the MCMC (use (3) 
for each lattice site x 

for each rectangular path r x t starting in x 
compute P rx t (x), the product of links along the path 
compute v = -log(P rx \x))/(t x L 4 ) 
add v to V(r) 
return V(r) 



Fig. 7. Example of algorithm to compute the static quark-antiquark potential. Notice the role of 
steps 7-9 is to create the initial configuration, step 11 loops over the Markov chain, step 12 creates 
the next configuration in the chain (using Metropolis, Gibbs sampling or other algorithm), steps 
13-16 measure the operator, and step 17 averages the results. 

For Q = SU(n) and n > 2 there is no exact technique but a common recursive 
technique 21 consists of 

U(x,u)^U'(x,u) = l[G ij (41) 

where Gij are random SU(2) matrices that acts on the ij subgroup of SU(n). 
The general algorithm is reported in fig. 6. 

3.3. Example: Quark- Antiquark Potential 

According to Quantum Mechanics, a state consisting of one static quark and 
one static antiquark separated by a distance r evolves in time by aquiring a phase 
e lHt where H, the Hamiltonian, is equal to the static potential between the quarks 
H = V(r). In fact, kinetic contribution to H is zero because the quarks are static. 

On a Euclidean lattice the same state evolves according to e' v ^ t , therefore the 
potential V(r) can be determined by computing the lattice expectation value of the 
operator that corresponds to this system. 

The quark and the anti-quark have to be created at the same point, separated 
at distance r, evolve for a time t and then reunite and annihilate. Due to the fact 
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Fig. 8. The static quark potential for quenched and two flavor dynamical QCD for different 
masses of the sea quarks 22 . Results are in units of the lattice spacing. 

that an antiquark is nothing other that a quark propagating backwards in time, the 
operator that corresponds to this system is the product of links around a square 
path of size r x t 

p rxt (x) d ^ U(x,ti)U(x + jj,,fi)...U{x + (r-l)fc,n) x 

U[x + rjl, v)U{x + rjj, + £>, u)...U(x + rp, + (t — l)z>, v) X 

U(x + r(i + ti>, —fj,)U(x + (r — l)/t + ti>, —n)...U{x + /t + (t— l)i>, —fj) x 

U(x + tP, -v)U{x + (t- l)u, -u)...U{x + P, -v) (42) 

The static potential can therefore be determined by measuring the left hand side 

of 

(0|^P rxt (x)|0) oce- v(r) * (43) 

X 

that is measuring 

V{r) = - j log I i- > ] ReTrP r x 1 (x) I ( 44 ) 



f^E JteTrprxt ( a; )j 



Fig. 8 shows the result of this computation 22 . The open squares represent the 
static potantial for a pure SU(3) gauge field theory in 4D (s = 4) referred to as 
quenched QCD. Notice that, for long distance, V(r) ~ err where a is the string 
tension. The solid points in figure represent the same static potential with a gauge 
field coupled to two dynamical light quarks for two different values of the quark 
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mass, to. We expect a change of regime when V(x) > 2m. In fact, according to cur- 
rent models of confinement 23 , when the static potential exceeds the energy required 
to create a new couple of quark and antiquark, these new particles pop up from 
the vacuum and screen the interaction between the original static quarks. The plot 
shows an indication of this phenomenon. 

4. Fermions 

Fermions, here identified by tp, are solutions of the Dirac equation. The latter can 
be derived from the continuum Euclidean action 

SE rmC = J d 4 x^ a (x)(^ fj D^ + m8 a pW {x) (45) 

where = d^ l + igA fl (x) is the covariant derivative, a/3 are spin indices and ij are 
gauge indices. 

From now on the gauge indices ij will often be omitted and we will restrict to 
4D, s = 4, since fermions are only defined for even number of dimensions. 

Notice that because of the Wick rotation, gamma matrices have to be rotated 
too. The Euclidean gamma matrices are hermitian (7^ = 7^) and satisfy {7^, 7„} = 
2<5 M „. Two possible choices for Euclidean gamma matrices in four dimensions are 
listed in the appendix. 

In the absence of gauge interaction, = and the simplest symmetrical 
discretized derivative looks like 

WW =' 5 [#c + £)-#«;-£)] (46) 

The most local generalization of eq. (46) that preserves the gauge invariance of the 
action is 

D^(x) = f \ [U(x, M Mx + £) - U\x - /i, ^{x - A)] (47) 
With this definition for the derivative, eq. (45) becomes 

s mrac = J- Mx)Q afj (x, x')Mx) (48) 

and Q is called fermionic matrix 

Q^ ve (x,x') = mS XtX ,S a0 +J2^f}l \. U ( X >V)S X +H, X ' - U{x,-ii)5 x -^ x ,] (49) 

There is still a problem with this naive action. The quark propagator S is defined 
as the inverse of the fermionic matrix Q, i.e. 

S%{x,x') de J (0|T{^(x'),CW}|0) - (Q-X^x') (50) 
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In absence of gauge field (U (x, /j,) = 1) for a massless quark (m = 0), for a single 
Fourier component of the fermionic field (ip(x) — e lp ' lX '') the inverse propagator 
reads 

S- 1 =Q = m + i^^sm{ Ptl ) (51) 

fj- 

(P/j is in units of 1/a here). This propagator has 16 poles as opposed to the single 
pole at p = for the continuum propagator. The additional poles arise when the 
spatial components of p are equal to ir. 

The physical interpretation of the additional poles is that this naive discretiza- 
tion of the action describes 16 degenerate fermions as opposed to a single one. Dif- 
ferent solutions to this problem lead to different implementations of lattice fermions. 
We consider here Wilson, clover, staggered, overlap, and domain wall fermions. 



4.1. Wilson Fermions 

Wilson proposed to remove the additional poles by giving mass to the corresponding 
modes 19 . This is done by adding a new term to the Lagrangian proportional to 

r^ a {x)D^D^ a {x) (52) 

This corresponds to replacing cq. (49) with 

Q a p(x,x') = (m + Ar)8 x ^5 al3 - |X),J (r - ^) af3 U{x : n)5 x+jl , x , - 

{r + ^) af3 U{x,-ii)5 x -^] (53) 



Introducing the definition 



(54) 



2m + 8r 

and scaling the field ip, eq. (53) can be rewritten as (omitting spin indices) 
Q w (x, x') ^S xy -K^2 [(r - r)U{x, n)5 x+fi , x - ( r + -fMx, -/i)5 x - A ,^] (55) 

which is the standard form for the Wilson fermionic matrix. Notice that any value 
of r > will do the job and one usually chooses r = 1. The choice r — corresponds 
to the naive action eq. (49). 

In the Wilson fermionic action the fermion mass is traded in for the adimensional 
parameter k defined in the asymptotically free limit by eq. (54). In presence of 
gauge interaction k is renormalized and eq. (54) holds only approximatively. This 
renormalizion shifts the value of k that corresponds to chiral fermions (m = 0). From 
now on we will identify with k* that value of k that corresponds to a chiral fermion. 
k* is not known a priori but it can be determined numerically as it corresponds to 
a pole of Q^ 1 . 
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For Wilson fermions on a cold configurations, U = 1, the fermion propagator 
reads 

S- 1 = Q = to + iJ2 7 M sin(p M ) + 2r ^ sin 2 (p„/2) (56) 

The main problem with Wilson fermions is that the fcrmionic matrix for to = 
does not anticommute with 75 and therefore the action is not invariant under the 
global chiral symmetry 

il>{x) -» e l75 V(x) $(x) -» V5(x)e l7Be (57) 

which is a symmetry of the continuum Dirac action. 

The breaking of chiral symmetry invalidates chiral perturbation theory which 
is used to guide guide the extrapolation of the spectrum to the limit to — > 0. This 
extrapolation is crucial in order to compute correlations involving fermions with 
mass smaller than 1/a. 



4.2. Clover Fermions 

The simplest way to Symanzik improve Wilson fermions is by shifting fermions 
according to 

[l + crfD^x) (58) 

and tune c in order to cancel any O(a) dependence in the correlations. The effect 
on the action is the same as replacing 24,25 

Q sw (x,y) = Q w (x,y) ^^W^^ ( 59 ) 

where F^ v = [D^,D V ] is a lattice version of the electromagnetic tensor which, in 
terms of the links, can be expressed as 

F\xv = (B- Bt)/ 8 
B = P(x, n, v) + P(x, /i, -v) + P(x, -fi, v) + P(x, -n, -v) (60) 

Eq. (48) with (59) is referred to as Sheikoleslami-Wohlert (SW) action or simply 
clover action because of the expression for F, eq. (60). 

The value of the coefficient csw does not affect the results in the continuum 
limit a — > but does affect the rate of convergence in this limit and the symmetry 
restoration for those symmetries broken by the lattice. 

There are three standard techniques to choose the parameter csw 

• 1-loop improvement 26 



c sw = 1 + 1.5954//? 



(61) 
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• Tadpole improvement 27 (i.e. rcsumming the contribution of all tadpole 
graphs to the renormalization of the tree-level csw)- 

c sw = 4 (62) 

where Uq is the average of ^ReTrP(x, \i,v) and it can be extracted from 
numerical simulations. 

• Non-perturbative improvement 26 . This is the most sophisticated technique. 
The idea behind it is that of determining the improvement coefficients for 
the different operators (including csw) by measuring independently on lat- 
tice the left and right-hand side of Ward identities and imposing the con- 
straint that they coincide. The results for csw can be summarized by the 
following fitting function (valid only for [3 > 5.7) 

/3 3 - 3.648/3 2 - 7.254/3 + 6.642 
Csw = /33 - 5.245S/? 2 (63) 

Even if different techniques give different results for csw they are consistent with 
each other provided the operators are improved by adopting the same procedure. 
Whichever improvement technique is used, the SW action generates correlation 
amplitudes that converge to the continuum limit up to correction of the second order 
in a and order 1 in g 2 (for 1-loop) or exactly (for non-perturbative improvement). 



4.3. Heavy Fermions and Fermilab Action 

As mentioned previously the lattice regularization acts as infrared cut-off and pre- 
vent particles with mass mq higher then 1/a to propagate properly. This presents 
a problem for the simulation of heavy quarks since in typical LQCD computations 
the lattice spacing is of the order of (2GeV) _1 . There are four ways to implement 
heavy fermions on the lattice. 

Extrapolation: Perform simulations with fermion masses lighter than the cut- 
off, m < 1/a, and extrapolate at the physical heavy fermion mass m — > ffiQ. 

Static fermions: Perform simulations in the static limit, wiq — -> oo. In this 
limit the Wilson action becomes the HQET 30 action with fermionic matrix given 

by 

Q HQET {x, x') = l —^U{x, 0)5 x+6jX + 1 -^U(x, -0)S x _ 6jX , (64) 

In this limit the a fermion propagator is known exactly and it is the product of links 
in the time direction 

S(x, x') = —^U{x, 0)U(x + 0, 0)...U(x' - 0, 0)S x ^6(x - x' ) + 

^L(U(x, 0)U(x + 6, 0)...U(x' - 6, 0)?6 x ^6(x' - .T ) (65) 
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Notice how static fermions require only one spin component because no term in the 
action mixes different spin components. Computations in the static limit are also 
useful to guide the extrapolation in the approach previously mentioned. 

Non relativistic limit: Adopt a non-relativistic approach (NRQCD) and per- 
form an expansion of the fcrmionic action around the on-shell momentum 31 of prop- 
agating fermions — p° n - shc11 + k^. Non-relativistically p° n " shc11 = ttiqv^ where 
is the on-shell velocity of the heavy fermion and is the off-shell momentum. 
In QCD, fc M is of the order Aqcd- This expansion results in the NRQCD action 
described by the fermionic matrix 

Q NRQCD (x, x') = i 7 (D + v D + 0(l/m Q )) (66) 

NRQCD fermions require two spin components which are mixed by O^/ttiq) terms. 
Corrections can be taken into account systematically. 

Fermilab action: The Fermilab action is an effective action that interpolates 
smoothly between the regular fermions and the static limit, and eliminates errors 
proportional to (arriQ) n . This is achieved by taking Wilson fermions with the clover 
action and introducing different couplings for space-like and time-like interaction 
terms in the Lagrangian. Moreover, in the spirit of Symanzik, higher order operators 
are added to the Lagrangian to cancel discretization terms that become sizable for 
large masses. A Discussion of these operators up to dimension 4 in Aqcd/itiq 
(HQET) and dimension 8 in (NRQCD) can be found in 28 - 29 . 

4.4. Staggered Fermions 

The Kogut-Susskind fermions, also known as staggered fermions provide an alter- 
native to Wilson fermions. 

The idea of Kogut and Susskind is that of interpreting the 16 poles of the naive 
fermion propagator, eq. (51), as due to the 4 spin components of 4 different types 
of fermions, here referred as flavors, which live on a blocked lattice 32 ' 33 ' 34 ' 35 ' 36 as 
shown in fig. 9. 

In order to introduce staggered fermions we will adopt the following notation: 

• x, x' will indicate points on the full lattice. 

• y, y' will indicate points on the blocked lattice 

• z will label the vertices of a hypercube of side a so that each x has a unique 
representation as y + z. <E 0, 1 are the four-spacetime components of z. 

• ipaa{y) will represent the four fermion flavors where a is the flavor index 
and a is the spin index. Note that ip is now defined on the blocked lattice. 

• \{x) will indicate the proper staggered field which corresponds to ip but is 
defined on the full lattice. 

• We will omit the color index since the formalism in completely transparent 
to it. 
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The map between naive lattice fermions and staggered fermions is realized by 

= n (y + * W(y> v + z )x(y + (67) 

z 

where 

^ {x) d g (68) 

and Vd/j-y + z) is a product of links connecting y to y + z that makes ip(y) transform 
correcly under local gauge transformations. The map in eq. (67) between ip <-> x is 
invertible. 

The naive action of eqs. (48,49) after substituting in eq. (67) becomes 

S E = £ <MlWA. + ™)<Ml/) 
y 

= ^{y + ^iy + z^i-i^D^ + m^iy + z^xiy + z') 
y,z,z' 

= E !Z ^[x t (zM*, M )xOr + /i) - X^M*, - M )x(z - /i)] + 
^m|x(a;)| 2 

^^(ilQ^t^lx^') (69) 

where 

r?(x, M ) = Itr {fit^h^C* ± M )} = (_i)E„<„ (70) 

and 

Q KS, (x,a;') = m5 XiX / + ^y^vjx, n) [U(x, m)<Wi,x' - U(x, — /u)5 x -i,x'] (71) 

Notice that in deriving the staggered action one uses a derivative term 
defined on the full lattice and not on the blocked lattice. This has the effect 
of coupling the different fermions in nontrivial ways and breaks the continuum 
SU {A) fi avor x SU{4) sp i n symmetry down to the discrete subgroup SWa x T4 where 
SW4 is the hypercubic subgroup of Euclidean rotations, SO (A), and is the Clif- 
ford algebra generated by the 7 matrices. This symmetry breaking allows us to 
identify the different flavors. 

In fact, in the naive discretization of the action (45) the doubling problem is 
related to the lattice symmetry 37 

fa) fa) = ]\{i l5l ^e^fa) (72) 

where z is in the hypercube. Eq. (72) suggests that the naively discretized eq. (45) 
contains interaction terms that couple one fermion mode into another by emission 
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of a hard gluon. For staggered fermions these modes correspond to different flavors 
and, therefore, can be distinguished. 

While Wilson fermions completely break the axial symmetry, staggered fermions 
break SU{4)v x SU{&)a only partially and the subgroup U{l)v x U(1)a remains 
unbroken thus preserving some form of the PCAC relations 38 which are important 
to guide the extrapolation to the chiral limit of quantities of phenomenological 
interest involving light fermions. 

Another advantage of staggered fermions over Wilson fermions is that the inver- 
sion of the corresponding fermionic matrix Q, i.e. the computation of the fermion 
propagator, is about eight times faster for the same lattice size. 

One disadvantage of staggered fermions is that they entangle spacetime and 
flavor symmetries thus making it difficult to engineer physical states of definite 
quantum numbers. 

In terms of the spinors the most general pseudoscalar meson can be written as 

^(y) = e b ^a, a (y) T ^Mv) (73) 

where T is 7 s or 7°7 5 . 

For this operator to have the quantum numbers of a pion, £ ab must be an element 
of the 15 representation of SU(4)fi avor . The choice £ ab = 5 a b would correspond to 
the singlet, the 771 . 

Different choices of basis for the £ oh matrices are present in the literature and 
they are equivalent to each other up to an unitary transformation since there is only 
one irreducible representation of the Clifford algebra in dimension 4, the algebra of 
ordinary gamma matrices (associated to the group T^). Therefore, according with 
usual conventions, we choose the following basis for the £ matrices 

£(5) = (7 5 )*, £(„) = (YT, = (tV)*, £(„,) = (7"7T (74) 




Fig. 9. Graphical representation of a 2D lattice and its blocked lattice (left). The center figure 
represents an example of an incorrect 2-point correlation amplitude between staggered fermions. 
The right figures represents two correct 2-point correlation amplitudes corresponding to different 
blockings of the same lattice. 



Fig. 9 represent a lattice, a blocked lattice, and how staggered mesons may 
propagate from one hypercube to another hypercube. Notice that source and desti- 
nation hypercube has to be consistent with the same blocking of the lattice or the 
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meson propagator will correspond to a nontrivial mix of the different mesons and 
will exhibit an oscillating behavior. 

Staggered fermions with the action eq. (69) converge quadratically to the con- 
tinuum. One way to Symanzik improve this action and achieve 0(a 2 ) improvement 
is by replacing each link that appears in eq. (71) with a weighted sum of links and 
paths connecting the same endpoints as the link. The choice of the weigth factors 
that cancels all 0(a 2 ) effects is called ASQTAD action 39 . 



4.5. Chiral Fermions: Overlap 

In continuum space the massless fermionic matrix is Q = 7^(9^ + igA^) and it 
satisfies the chirality condition 

Ql 5 + 7 5 Q = (75) 

therefore eigenstates of Pl = (1 — 7 5 )/2 and Pr = (1 + 7 5 )/2 are not mixed by 
the action. The only term in the action that couples L and R chirality states is the 
mass term. 

Eq. (75) docs not hold for Q w , Q sw and Q KS and, in fact, the Nielsen-Ninomiya 
theorem 40 states that it is not possible to define a local, translationally invariant, 
hermitian Q that preserves chiral symmetry and does not cause doublers. Never- 
theless one may look for a fermion formulation on the lattice that exhibit chiral 
symmetry and no doublers by requiring that chiral symmetry holds for on-shell 
states only. For on shell states chiral transformation can be written as 

^ ^ e tfV(l-aQ/2ty ^ ^ ^7 5 d-aQ/2) (?6) 

In fact on shell Q-ip — and eq. (76) becomes the prdinary chiral symmetry tran- 
sormation. Requiring that the fermionic matrix Q be invariant under infinitesimal 
transformation leads to the Ginsparg- Wilson relation 

Q 7 5 + 7 5 Q = aQj 5 Q (77) 

Neuberger 41 ' 42 proposed the first fermionic matrix Q that satisfies the Ginsparg- 
Wilson relation and thus provides exact chiral symmetry on the lattice 

QOverlap = 1 + y^g^S (Q^ _ ^ ^ 

This formulation is referred to as overlap fermions. 

Since any hermitian matrix X can be written as £ = ADA^ where D is a diagonal 
matrix, one can define any function / of a hermitian matrix via 

/(E) = A/(D)At (79) 

and f(D) is a diagonal matrix with diagonal elements given by f(D)u = f(Du). 
The argument of sign function in the definition of the overlap fermionic matrix is 
hermitian therefore eq. (78) is well defined. 
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The sign function can be approximated by polynomials and, in general, it is more 
difficult to deal with than Wilson, clover or even staggered fermions. All numeri- 
cal techniques 47 for calculating a det Q° verla P a s required for dynamical fermions 
and for calculating (QOveriap^-i can ^ e v j ewec j as different ways of approximating 
cq. (78). 

4.6. Chiral Femions: Domain Wall 

An alternative way to implement lattice chiral fermions was developed by Kaplan 43 
in 1992. It is referred to as using Domain Wall (DW) fermions 44 ' 45,46,47 . This for- 
mulation consists of replacing each fcrmion ip a (x) with multiple fermions labeled 
by a new index k 

i> a (x) -» * a ,k(aO (80) 

The DW action is 

sE w = T,*aAx)Q2%kk>M*f>,k>&) (si) 

x.k 

where 

Q™ k ,(x,x')=Q™(x,x')5 kk ,+ 

(P L ) af 30(k <N 5 - l)5 k +i, k < + {P R ) a pO{k > 0)5 k -i,k> + 

(PL)af3$k,N 5 -l<>k' ,0 + {P^a^kfi^k' ,JV 5 -1 (82) 

&5 and N§ are parameters of the model. The index k runs from to N$ — 1 and it 
is usually interpreted as a 5th dimension of the system. Notice that the gauge field 
that appears in Q w is independent on this 5th dimension. 

The computation of ip' = Q~ 1- ip for a regular 4 dimensional input fermionic field 
ip is performed in three steps: 

• The input field ip is mapped into the 5-dimcnsional DW field \&. L compo- 
nents are mapped into the k — slice and R components are mapped into 
the k = N 5 — 1 slice. 

y a ,k( x ) =° for < k < N 5 - 1 

*a,JV 5 -l = (PR) a p^p(x) (83) 

(84) 

• The DW fermionic matrix, eq. (82), is inverted numerically using one of the 
algorithms explained later 

*' = (Q^r 1 * (85) 

• The output DW field is mapped back into the output 4 dimensional field 
V'by 

1>' a {x) = (P L U% fi (x) + (PrU^n^x) (86) 



29 



The effect of the DW action is that of projecting the L modes of the 4D fermion 
to one wall and the R modes to the other wall, thus different chiralities are mixed 
only by the explicit mass term in Q DW , eq. (82). The use of the Wilson action in 
each slice of the extra dimension k guarantees that no doublers are present. 

It is possible to translate the DW fcrmions into an equivalent four-dimensional 
approximation for the Neuberger operator 48 , eq. (78) and therefore Q° verla P \ s 
a local operator. In fact, any rational polynomial approximation of the overlap 
operator is equivalent to a domain wall formulation with a finite 5th dimension 49 . 



4.7. Quenched and Dynamical Fermions 

Including fcrmionic field variables in the MCMC configurations is not practical. 
The usual technique to deal with fermions is integrate them out analytically so that 
fermions do not appear in the PI measure 

Let's consider a system consisting of a gauge field coupled to n/ degenerate 
fcrmions 

S E = S 9 E au9e + n f S E erml 

= 2^ E ReTrP{x,u,u)+n f J2^(x)Q a p{x,y)ij {y) (87) 

x,u=£v x.y 

For a typical correlation amplitude, fermions can be integrated out as follows 

(0| ...M X )M X ')- |0) - J [d[/][dV]...VaW^(^)-e- 5r " e - n/Srmi (88) 
= J [dC/]...Q^(x,a;')...e^ a " 9e+ ^ losdctQ (89) 

where the field configurations are now generated at random by sampling from a 
probabily distribution proportional to cxp(—S E uU ) with 

s full = ggauge _ ^ {Qg ^ Q ^ 

In case the correlation amplitudes involve multiple possible Wick contractions, 
one has to sum over all possible Wick contractions. 

(o| ...M*)M*%-MJ')M*'") l°> = ^ E -Qap( x > x ')-Qi5( x " > x "')- + 

u 

...Q-l(x,x'")...Q-^x",x')... (92) 

For staggered fermions, since Q KS represents four degenerate flavors as opposed 
to a single one, therefore eq. (91) must be replaced by 

s fullKS = ggauge _ ^ ^ qKS ^ 
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It is still debated whether the above procedure is correct, since there is no 
known local operator that corresponds to the 4th root of Q KS . Anyway, there are 
indications from numerical studies 50 in 2D that the following relation may hold in 
4D 

(det Q KS ) 1/4 = det Q ^ rla P + 0{a 2 ) (94) 

which would provide a solid theorical justification for eq. (93). 

From the definition it is evident that Q is a sparse matrix. For Wilson fermions 
its dimensions are M x M and M = 8n(L/a) 4 (real and imaginary part x 4 spin 
components x n color components x number of lattice sites); Q has diagonal ele- 
ments set to 1 and has only 4 other elements different from zero for each row and 
column (corresponding to +/z and —fi). For staggered fermions Q KS is 16 times 
smaller because it has no spin indices and this makes its inversion much faster. For 
domain wall fermions Q DW is iVf larger than Q w thus making the inversion of 
domain wall fermions much slower than for Wilson fermions. 

One approximation that has been used and abused consists of setting rif = in 
the full action. This approximation is called quenching. It has the effect of ignoring 
second quantization for fermions. This is equivalent to, in a perturbative language, 
ignoring fcrmion loops. Quenching introduces unknown systematic errors in the 
computation and its only justification is that the computation of det Q is non-local 
therefore generating the Markov Chain with n/ ^ is very computing intensive. 

For some quantities such as the static quark potential, fig. 8, the effect of quench- 
ing is very small. For other quantities such as the light spectrum of QCD, its effect 
is sizable, fig. 16. Quenching also affects the behavior of the spectrum in the chiral 
limit as explained in 51 . 

The contributions det Q to the action is referred to as dynamical fermions or 
sea quarks in the context of LQCD. 

4.8. CPTH Theorem on the Lattice 

Lattice correlation amplitudes computed using the action in eq. (87) are invariant 
under charge conjugation, C, parity, P, time reversal, T, and a new symmetry, H. 
These symmetries apply to the measurement of an operator on each gauge config- 
uration U. 

It is useful to write how these symmetries affect a fermion propagator S = Q^ 1 
and make explicit the dependence on all indices and on the gauge configuration U. 

• Charge conjugation, C: 

S%(x,y,U) = (7VW^V(2/,^ C )(7 2 7°W (95) 

• Parity, P: 

S%(x,y,U)= 1 l a ,S i J, ,(x p ,y p ,U P h% (96) 
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Algorithm: Minimal Residual Inverter 
Input: ip, Q 
Output: V' 

Temporary fields: q,r 

r = tp — Qtp 
tp' = tp 
do 

q = Qr 

a=(q- r)/(q ■ q) 
tp' = tp' + a r 
r = r — a q 
residue = r ■ r 
while residue > precision 



Fig. 10. Minimal Residual, a numerical algorithm to compute ip' = Q 1 tp by minimizing Qip 1 = 

• Time reversal, T: 

S%(x,y,U) = (7VW^,(x T ,y T ,^ T )( 7 VW (97) 

• H symmetry: 

S%{x,y,U) =ll al SZ p ,{y,x,U) 1 % a (98) 

U p , U c , U T are the parity reversed, charge conjugate, and time reversed gauge 
configurations respectively. 

4.9. Inversion Algorithms 

One way to invert the matrix Q is by using a stochastic technique. 
In fact for any hermitian positive definite matrix £ the following exact relation 
holds: 

Sy 1 = \ J d^o...d^_i0;^e-*» E »™*™ (99) 

And substituting in the above equation S = Q^Q and multiplying both terms by 
one obtains 

(Q-^t^x 1 ) = iy M(Q<t>%'tf)4 l (x)e-+W*» (100) 

This multidimensional integral can be computed via Monte Carlo. The field <p in 
this context is usually referred to as pseudo-fermionic field. The stochastic com- 
putation of the full inverse propagator must be done for each gauge configuration 
and therefore it is computationally expensive. Various techniques for reducing the 
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1 Algorithm: Stabilized Biconjugate Gradient 

2 Input: ip, Q 

3 Output: ip' 

4 Temporary fields: p,q,r,s,t 

5 

6 Ip' = 1p 

7 r = ip — Qip 

8 q = r 

9 p = (zero field) 

10 s = (zero field) 

11 p = p' = a = cj = l 

12 do 

13 p = q ■ r 

14 j 8=(p/p')(a/w) 

15 p' = p 

16 p = r + [3p — (3ljs 

17 s = <5p 

is a = p/(g • s) 

19 r = r — a s 

20 t = Qr 

21 uj = (t ■ r)/(t ■ t) 

22 ip' — ip' + uj r + a p 

23 residue — r ■ r 

24 while residue > precision 



Fig. 11. Stabilized Biconjugate Gradient, another numerical algorithm to compute ip' = Q 1 ip 
by minimizing Qtp' = ip. 



variance in the integration and reduce the statistical noise have been proposed by 
many authors 52 ' 53 . 

In most cases one does not need the full inverse Q^ 1 and it suffices to solve in 
ip' the equation 

Qtp' = ip ip' = (101) 

for a small set of given input vectors ips. This inversion can be performed by mini- 
mizing the norm of the residual vector 

r = ip-Qip' (102) 

The two algorithms commonly used for performing this numerical minimization are 



33 



the Minimal Residual, fig. 10 and Stabilized Biconjugate Gradient fig. ll c 

4.10. Dynamical Fermions Algorithms 

For any Hermitian matrix S 



Thus for two degenerate flavors the contribution to the action due to fermions is 



where <p are pseudobosonic fields. Eq. (105) can be evaulated numerically using 
Monte Carlo. 

Various techniques based on the above equation have been developed to speed 
up the MCMC and to take into account odd numbers of dynamical quarks. 

The most common MCMC algorithm for dynamical fermions is the Hybrid 
Monte Carlo algorithm discussed in 54 > 55 . 

The technique used in 22 consists of observing that dct Q^Q = (det j 5 Q) 2 where 
j 5 Q is hermitian, and approximating the determinant with a truncated determinant, 
defined as the product of the eigenvalues of 7 5 Q below some cut-off. The eigenvalues 
arc computed by diagonalizing 7 5 Q using the Lanczos algorithm. 

One of the most promising techniques in terms of efficiency for light fcrmion 
mass is a version of the Hybrid Monte Carlo based on the Schwartz Alternating 
Procedure discussed in 56 and 57 . 

Theoretical work on algorithms for dynamical overlap fermions are proposed in 
58 and 59 . Some preliminary phenomenological results can be found in 60 . Algorithms 
for dynamical domain wall fermions are discussed in 61 and 62 . 

4.11. Example: Pion Mass and Decay Constant 

We consider here, as an example, the computation of the mass of the pion, , 
and the pion decay constant, /„., for a SU{n) gauge theory (for n = 3 this is QCD) 
as defined by the action in eq. (91). 

In order to be able to put a pion on the lattice one needs to have a definition 
of the former, i.e. one needs to have an operator function of the gauge and quark 
fields that has the pion as lowest eigenstate. 

c In writing fermionic algorithms we adopted the following notation: 




(104) 



(detg) 2 =dctQt Q = J_ / [dfle-WQr 1 * 



(105) 



E tf*w = cw 



2 



E cw=cw 



(103) 



y,0,j 



34 



Algorithm: Compute the static quark-antiquark potential 
Input: /3, k, size of gauge group n, number of MCMC steps N 
Output: C*(t) 

Create local array C n (t) and initialize it to zero 

for each lattice site x 
for each direction ^ 
set U(x, ^) to a random element of the gauge group SU(n) 

for c = 1 to N 

replace U with the next configuration in the MCMC 
for each spin component a 
for each color component i 
make a field ^(x) = 
for each lattice site x on timeslice x° = 

set ip a - i (x = 0) = 1 
compute ip' = Q _1 V> 
for each lattice site x 

add \iP'(x)\ 2 \o C 7r (t = x°) 
return C n (t) 



Fig. 12. Example of algorithm to compute a pion progator. Notice the role of steps 7-9 is to 
create the initial configuration, step 11 loops over the Markov chain, step 12 creates the next 
configuration in the chain (using Metropolis, Gibbs sampling or other algorithm), steps 14-18 
measure the operator, and step 19-20 average the results separately for each lattice time-slice. 



The pion is the lightest pseudo-scalar in the theory and the simplest operator 
that transforms as a pseudo-scalar under rotations is 



*.(x)^]T^*) Q75 Q %0r) (106) 



We identify with \Ek) its eigenstates and Ek the corresponding ordered eigenvalues 
so that, by definition, E is and \E ) is \n). A quantum mechanical computation 



35 



0.22 — i — i — i — i — | — i — i — i — i — | — i — i — i — r 




0.00 0.05 0.10 0.15 

(m,+m y )r 1 x (Z m /ZT) 

Fig. 13. Global fit of partially quenched data from the MILC collaboration 63 and chiral extrap- 
olation of /V for various input values of the propagating and dynamical quarks' masses. Here m x 
and m y represents the bare u and d quark masses respectively. 



shows that 

def 



CAyo - x ) = J FT°FT° (0| *(x)^(y) |0) (107) 
= ]T FT® FT® (0| V?(x) \E k ) -i- (E k \ ¥(y) |0) 



2E k 



|2 



2E k 

= l(oi^(Q)K)l 2 c _, ro , fan _ zn) | ^ 

_ Jir m ir c -i m ,( w -m) _|_ ^ (108) 

Here FT® — ^ x is the Fourier transform in the spatial components of x at 
zero momentum. The dots indicate exponential terms that oscillate faster that the 
leading term. The last step follows from the definition of f w , the pion decay constant. 

After a Wick rotation 

CV(y - x Q ) = £dpL e -iMi«-xo) + ... (109) 

and the fast oscillating terms, represented by the dots, are replaced by fast decaying 
exponentials. For tjq — xq = t and considering periodic boundary conditions 

C„(t) = tiHhL[ e -™At) +e -m„(L-t)] + ... (n0) 
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Fig. 14. Example of currents used on lattice and their relative quantum numbers, tp, -0' and -0' 
are different flavors. The superscripts i, j and k are gauge indices. 



The lattice formulation of QCD provides the numerical technique to evaluate 
the left hand side of eq. (108) 

C„(t) = FT^ J[d(f ) }[dA}^(x)^(y)e- SE 

- Jf E FT°FT°ReTr^Q-\x, y)^Q-\y, x)) 
u 

= ^EEEi^" 1 ^y)i 2 ^-iyo-^oi) an) 

(7 £C y 



In the last step we used H symmetry on the second propagator. 

By computing eq. (Ill) for different values of t and fitting the results with 
cq. (110) one can extract both m v and f^. Some numerical results 64 for / w computed 
from eq. (Ill) for different values of the quark masses are shown in fig. 15. The 
extrapolated f n is compared with the experimental results. 

In practice, because of dimensional transmutation, one always computes pure 
numbers such as m„ in units oil/ a and one can eliminate such dependence on a 
by computing ratios of masses or other dimensionless ratios. 

Similarly one can compute masses and decay constants of other particles by 
choosing the right operator. A list of operators for various interesting states is 
listed in fig. 14. For a deeper analysis on how to built this type of operators for 
baryons can be found in 65 . 

Fig. 15 shows the computation of the neutron mass for different fcrmion formu- 
lations at different lattice spacing 66 . As expected, within error, they all agree with 
each other in the contiuum limit. 
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quenched comparison 
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0.1 0.2 0.3 0.4 0.5 
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Fig. 15. The plot shows the extrapolation to the continuum limit of the neutron mass computed 
in LQCD using various types of fermions 63 . Improved glue refers to an 0(a 2 ) improved gauge 
action. Different formulations agree with each other within the error but clover and staggered 
converge faster then Wilson as expected, mps in the plot refers to the pion mass which is used to 
set the unit scale r\. 



5. Error Analysis 

There are two types of errors in any LQCD computation. Statistical errors and 
systematic errors. Statistical errors are under control today. The main source of 
systematic error are discussed below, they are being addressed by recent computa- 
tions, and will be removed in the near future. 

• Discretization. The typical lattice spacing is today of the order of (2GeV) _1 . 
This introduces discretization errors that are sometime difficult to quantify. 
One effect, for example, is the breaking of continuous rotational symmetry. 
One way to reduce discretization effects is by means of Symazik improved 
actions. 

• Quenching. This error is the most difficult to quantify and it is now being 
eliminated thanks to new algorithms and cheaper computing power. A great 
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Fig. 16. Recent LQCD results from the MILC collaboration 68 . The two images show ratios of 
LQCD results over experimental results for various QCD quantities of phcnomcnological interest 
without and with dynamical fermions, respectively. Dynamical fermions are implemented using 
0(a 2 ) improved staggered fermions. 



deal of effort has been put in this direction in the recent years. Recent 
computations with dynamical quarks have been a success but the dynamical 
quark masses are still larger then the physical u and d quark mass. 

• Chirality. Both Wilson and Staggered fermions break chiral symmentry and 
do not allow computations at zero quark mass. Moreover, the more the chi- 
ral limit is approached, the more expensive it is to invert the fermionic ma- 
trix Q. Alternative fermionic discretizations such as domain wall fermions 
and overlap fermions promise restoration of the chiral simmetry in the con- 
tinuum limit and provide a better approximation of continuum physics. 

• Matching. Most experimental quantities are usually expressed in the MS 
scheme, while the lattice computations are performed in a different regular- 
ization/renormalization scheme. The procedure to relate one to the other 
is called matching and it mainly perturbative in nature. When lattice re- 
sults are converted in the MS scheme using one-loop matching they become 
affected by matching errors as big as 20%. The computation of matching 
beyond one loop is a very challenging task. Some attempts to automate 
this process via a numerical approach to lattice pertubation theory can be 
found in 67 and 68 . It is important to stress that the use of the MS scheme 
is a convention and this matching would not be necessary if the lattice 
regularization were used everywhere. 

In LQCD, errors due to contiuum extrapolation, extrapolation to physical 
masses (both for heavy quarks and light quarks) and finite lattice size are today 
very much under control and below 2% for most quantities of phenomenological 
interest. 
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Fig. 17. Unquenched hadron masses and splittings compared with experimental values 69 . The T 
and cc columns are differences from the ground state masse as computed by HPQCD and Fermilab 
groups. The it, K and T(IP-IS) masses are used to fix the quark masses and the lattice spacing. 



Fig. 16 and fig. 17 show the effect of quenching on the light quark spectrum and 
how the latest dynamical LQCD computations 69,70 agree with experiment. 



6. Conclusion 

The lattice formulation provides a way to regularize, define and compute the Path 
Integral in a Quantum Field Theory. This formulation and the associated numerical 
techniques have been of crucial importance in deepening our knowledge of quantum 
field theories in physical regimes where perturbation theory fails, as in the case of 
QCD at strong coupling. For example, LQCD computations have played and are still 
playing an important role in understanding the process of quark confinement ' , 
determining the phase structure of gauge theories 73 , confronting QCD predictions 
with experiments 69 , and extracting fundamental physical parameters such as quark 
masses 74,75 and CKM matrix elements 76 ' 77 from experimental results. 

Three advantages of LQCD over most analytical techniques is that it is easier 
to understand, the effect of its approximations can be controlled numerically, and 
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the precision of any computation can be reduced arbitrarily just by increasing the 
dedicated computing time. 

Today most LQCD computations are limited by available computing power. 
In order to overcome this limitation many groups have been working on design- 
ing dedicated machines for LQCD such as APEnext 78 , QCDOC 79 and the Earth 
Simulator 80 . Other groups have focused on the development of software and al- 
gorithms optimized for commodity hardware such as PC clusters 81 and building 
infrastructures for the exchange of field configurations 82 . 

Most of the algorithms and the examples discussed here and many more are 
implemented in a free software library called FermiQCD 83 ' 84 ' 85 . Despite its name, 
FermiQCD is not LQCD specific but it is suitable for generic LQFT computations. 
It has a modular design and an easy to use syntax very similar to the one we have 
used in this paper (in fact, algorithms such as the MinRcs and the BiCGStab map 
line by line) . Moreover, all of the FermiQCD algorithms are parallel and they can 
run on distributed memory architectures such as PC clusters. FermiQCD has been 
used for production grade computations at Fermilab and other institutions. 

FermiQCD can be downloaded from www.fermiqcd.net. 
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Appendix A. Useful Formulas in 4D Euclidean Space-Time 
Appendix A.l. Wick rotation 

The Euclidean action is obtained from the Minkowskian one by performing a Wick 
rotation. Under this rotation the basic vectors of the theory transform according 
with the following table (E for Euclidean, M for Minkowski) 



E 


M 


E 


M 


x° 


ix° 


x l 


x l 


go 


-id 


Ql 


di 


A 4 


-iA 


A 1 


A l 


F°i 






Fij 


7° 


7° 




-irf 


7 5 


7 5 







and the integration measure transforms as follow 

expt-Sg) = exp(iS M ) (A.2) 
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where 



S E = j d 4 x E £. E [...} =-iJ d 4 x M jC M [...] (A.3) 

The choice d 4 x E = id 4 XM can be made, hence £e[...] — —£m[-~] 
The Euclidean metric tensor is defined as 

gg = -St" = diag(-l, -1, -1, -1) (A.4) 



Appendix A. 2. Spin matrices 

• Dirac matrices (Dirac representation) 

• Dirac matrices (Chiral representation) 



All the Euclidean Dirac matrices are hermitian. The following relations hold 

= 1 -{^,Y} = 5^ (A.7) 

^ = ^,71 (A.8) 

7 5 = 7°7 1 7 2 7 3 (A.9) 

and all the a pu are hermitian. 
Projectors 



• Traces 



1 - T 5 1 + T 5 



tr(YV) = 4<5^ (A.ll) 

trfr'W) = (A.12) 

trfr'WY 7 ) = 4(S^6 pa - 5 w 8 va + 5^5 pv ) (A.13) 

tr( 7 V7V7 ff ) = 4e^ pCT (A. 14) 



where e° 123 = -1. 



References 

1. M. Peskin and D.V. Shroeder, Quantum Field Theory, Addison Wesley (1994) 

2. C. Itzykson and J.-Z. Zuber, Quantum Field Theory, Mc Graw-Hill (1985) 

3. M. Creutz, Quarks, Gluons and Lattices, Cambridge University Press (1983) 

4. Rothe, Lattice Gauge Theories, World Scientific 



42 



5. I. Montvay and G. Miinster, Quantum Fields on a Lattice, Cambridge University 
Press (1994) 

6. R. Gupta, arXiv:hep-lat/9807028. 

7. T. Muta, Foundations of Quantum Chromo Dynamics, World Scientific (1987) 

8. W. Marciano and H Pagels, Phys. Rep. 36 (1978) 135 

9. L. Maiani and M. Testa, Phys.Lett. B245 (1990) 585-590 

10. S.P. Meyn and R.L. Tweedie. Markov Chains and Stochastic Stability. Springer- 
Verlag London, 1993 

11. G. Banhot, Rep. Prog. Phys. 51 (1988) 429 

12. J. Shao and D. Tu, The Jackknife and Bootstrap, Springer Verlag (1995) 

13. K. Symanzik, Nucl. Phys. B226 (1983) 187,205 

14. D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30 (1973) 1343 
D. J. Gross and F. Wilczek, Phys. Rev. D8 (1973) 3633 

H. Politzer, Phys. Rev. Lett. 30 (1973) 1346 
H. Politzer, Phys. Rep. 14C (1974) 129 

15. J. M. Rabin, Introduction to Quantum Field Theory for Mathematicians, IAS/Park 
City Mathematics Series, vol. 1 (1995) 

16. J. Glimm and R. Jaffe, Quantum Physics (second edition), Springer- Verlag (1987). 

17. M. Gockeler, R. Horsley, A. C. Irving, D. Pleiter, P. E. L. Rakow, G. Schierholz and 
H. Stuben, arXiv:hep-ph/0502212. 

18. Y. Aharonov and D. Bohm, Phys. Rev. (Ser. 2) 115 485-491 (1959) 

19. K. G. Wilson, Phys. rev. D10 (1974) 2445 

20. M. Creutz, Phys. Rev. D21 2308 (1980) 

21. N. Cabibbo and E. Marinari, Phys. Lett. 119B (1982) 387 

22. A. Duncan, E. Eichten and H. Thacker, Phys. Rev. D 59, 014505 (1999) [arXiv:hep- 
lat/9806020]. 

23. G. 't Hooft, Phys. Rev. (1978) 1 

24. B. Sheikoleslami and R. Wolhert, Nucl. Phys. B259 (1985) 572 

25. M. Luscher, Advanced Lattice QCD, Talk given at Les Houches Summer School in 
Theoretical Physics, Session 68, (1998); hep-lat/9802029 

26. M. Luscher, S. Sint, R. Sommer and H. Wittig, Nucl. Phys. B491 (1997) 344. 

27. CP. Lepage and P.B. Mackenzie, Phys. Rev. D48 (1993) 225. 

28. A. X. El-Khadra, A. S. Kronfeld, P. B. Mackenzie, S. M. Ryan and J. N. Simone, 
Phys. Rev. D58 (1998) 014506 

29. M. B. Oktay, A. X. El-Khadra, A. S. Kronfeld and P. B. Mackenzie, Nucl. Phys. 
Proc. Suppl. 129, 349 (2004) [arXiv:hep-lat/0310016]. 

30. M. B. Wise, arXiv:hep-ph/9311212. 

31. A. Manohar, Phys.Rev. D56 (1997) 230-237 

32. J. Kogut and L. Susskind, Phys. Rev. Dll (1975) 395 

33. T. Banks, J. Kogut, and L. Susskind, Phys. Rev. D13 (1996) 1043 

34. L. Susskind, Phys. Rev. D16 (1977) 3031 

35. G. W. Kilcup and S. R. Sharpe, Nucl. Phys. B 283 (1987) 493. 

36. M. F. Golterman, Nucl. Phys. B 273 (1986) 663. 

37. G. P. Lepage, Phys. Rev. D 59, 074502 (1999) [arXiv:hep-lat/9809157]. 

38. W. Lee and S. R. Sharpe, Phys. Rev. D 60 (1999) 114503 [hep-lat/9905023]. 

39. C. Bernard et al. [MILC Collaboration], Phys. Rev. D 58 (1998) 014503 [hep- 
lat/9712010]. 

40. H. B. Nielsen and M. Ninomiya, Nucl. Phys. B185 (1981) 20; Erratum Nucl. Phys. 
B195 (1982) 541 

41. H. Neuberger, Phys. Rev. Lett. 81, 4060 (1998) [arXiv:hep-lat/9806025]. 



43 



42. R. Narayanan and H. Neuberger, Phys. Rev. Lett. 71 (193) 3251; Nucl. Phys. B412 
(1994) 574; Nucl. Phys. B443 (1995) 305 

43. D. B. Kaplan, Phys. Lett. B288 (1992) 342 

44. Y. Shamir, Nucl. Phys. B406 (1993) 90 

45. V. Furman and Y. Shamir, Nucl. Phys. B439 (1995) 54 

46. V. Furman and Y. Shamir, Nucl. Phys. B 439, 54 (1995) [arXiv:hep-lat/9405004]. 

47. T. Blum et al., Phys. Rev. D 69, 074502 (2004) [arXiv:hep-lat/0007038]. 

48. A. Borici, Nucl. Phys. Proc. Suppl. 83, 771 (2000) [arXiv:hep-lat/9909057]. 

49. A. D. Kennedy, Nucl. Phys. Proc. Suppl. 140 190 (2005) [arXiv:hep-lat/0409167] 

50. S. Durr and C. Hoelbling, Phys. Rev. D 71, 054501 (2005) [arXiv:hep-lat/0411022]. 

51. S. R. Sharpe and Y. Zhang, Phys. Rev. D 53, 5125 (1996) [arXiv:hep-lat/9510037]. 

52. C. Michael and J. Peisa [UKQCD Collaboration], Phys. Rev. D 58, 034506 (1998) 
[arXiv:hep-lat/9802015]. 

53. A. O'Cais, K. J. Juge, M. J. Peardon, S. M. Ryan and J. I. Skullerud [TrinLat 
Collaboration], arXiv:hep-lat/0409069. 

54. R. T. Scalettar, D. J. Scalapino, and R. L. Sugar, Phys Rev/ B34 7911 (1986) 

55. S. Duane, A. D. Kennedy, B. J. Pendelton, and D. Roweth, Phys. Lett. 195B 216 
(1987) 

56. M. Liischer, Comput. Phys. Comm. 156 209 (2004) 

57. M. Liischer, JHEP 05 052 (2003) 

58. H. Neuberger, Phys. Rev. D 60, 065006 (1999) [arXiv:hcp-lat/9901003]. 

59. T. DeGrand and S. Schaefer, Phys. Rev. D 71, 034507 (2005) [arXiv:hep- 
lat/0412005]. 

60. Z. Fodor, S. D. Katz and K. K. Szabo, arXiv:hep-lat/0409070. 

61. P. Chen et al., arXiv:hep-lat/9812011. 

62. Y. Aoki et al, arXiv:hep-lat/0411006. 

63. C. T. Sachrajda, Lattice Simulations and Effective Theories, Lectures presented at 
the Advanced School on Effective Theories, Almunecar (Spain) June 1995; hep- 
lat/960527 

64. C. Aubin et al. [MILC Collaboration], Phys. Rev. D 70, 114501 (2004) [arXiv:hep- 
lat/0407028]. 

65. S. Basak et al., arXiv:hep-lat/0508018. 

66. C. T. H. Davies, G. P. Lepage, F. Niedermayer and D. Toussaint, Nucl. Phys. Proc. 
Suppl. 140, 261 (2005) [arXiv:hep-lat/0409039]. 

67. F. Di Renzo and L. Scorzato, JHEP 0410, 073 (2004) [arXiv:hep-lat/0410010]. 

68. H. D. Trottier, Nucl. Phys. Proc. Suppl. 129, 142 (2004) [arXiv:hep-lat/0310044]. 

69. C. T. H. Davies et al. [HPQCD Collaboration], Phys. Rev. Lett. 92, 022001 (2004) 
[arXiv:hep-lat /0304004] . 

70. C. Aubin et al, Phys. Rev. D 70, 094505 (2004) [arXiv:hep-lat/0402030]. 

71. C. W. Bernard et al, Phys. Rev. D 62, 034503 (2000) [arXiv:hep-lat/0002028]. 

72. M. Engethardt, Nucl. Phys. B Proc. Suppl. 140 2005 

73. J. B. Kogut and M. A. Stephanov, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 
21, 1 (2004). 

74. M. Gockeler, R. Horsley, A. C. Irving, D. Pleiter, P. E. L. Rakow, G. Schierholz and 
H. Stuben [QCDSF Collaboration], arXiv:hep-ph/0409312. 

75. C. McNeile, C. Michael and G. Thompson [UKQCD Collaboration], Phys. Lett. B 
600, 77 (2004) [arXiv:hep-lat/0408025]. 

76. D. Becirevic, arXiv:hep-ph/0211340. 

77. M. Okamoto, arXiv:hep-ph/0505190. 

78. F. Bodin et al. [ApeNEXT Collaboration], Nucl. Phys. Proc. Suppl. 140, 176 (2005). 



44 



79. P. Boyle et al, Nucl. Phys. Proc. Suppl. 140, 169 (2005). 

80. S. Aoki et al, Nucl. Phys. Proc. Suppl. 129, 859 (2004) [arXiv:hep-lat/0310015]. 

81. D. Holmgren, P. B. Mackenzie, D. Petravick, R. Rechenmacher and J. Simone, Pre- 
pared for CHEP'01: Computing in High-Energy Physics and Nuclear, Beijing, China, 
3-7 Sep 2001 

82. B. Joo and W. Watson, Nucl. Phys. Proc. Suppl. 140, 209 (2005) [arXiv:hep- 
lat/0409165]. 

83. M. Di Pierro, arXiv:hep-lat/0011083. 

84. M. Di Pierro, Nucl. Phys. Proc. Suppl. 106, 1034 (2002) [arXiv:hep-lat/0110116]. 

85. M. Di Pierro et al. [FermiQCD Colaboration] , Nucl. Phys. Proc. Suppl. 129, 832 
(2004) [arXiv:hep-lat/0311027]. 



